function printout_dispersion_eq(m, CONSTS, span, plot_data)

    omega_p = CONSTS.omega_p;
    omega_H = CONSTS.omega_H;
    num_points = 10; %30;

    omega_raw_vec = logspace(log10(0.01*omega_H)*(1+1e-5), log10(0.002*omega_p)*(1-1e-5), num_points)';
    omega_refined_vec = []; %linspace(1.24e+10, 1.28e+10, 30)';
    omega_vec = sort([omega_raw_vec; omega_refined_vec]);
    num_roots_vec = zeros(size(omega_vec));
    P_vec = zeros(size(omega_vec), 1);
    P_c_vec = zeros(size(omega_vec), 1);

    if (plot_data)
        figure;
        hold;
        title('p_{m,n} from {\omega}/{\omega_H}')
        xlabel('');
        ylabel('');
        %ylim(span);
        h_dispers = gca; 
        legend_vec = {};
    end

    for i = 1:size(omega_vec, 1)
        fprintf('i =  %f \n', i);
        fprintf('omega =  %f \n', omega_vec(i));
        LOCAL_CONSTS = CONSTS;
        LOCAL_CONSTS.f = omega_vec(i)/(2*pi);
        LOCAL_CONSTS = update_CONSTS(LOCAL_CONSTS);

        P_vec(i) = LOCAL_CONSTS.P;
        P_c_vec(i) = LOCAL_CONSTS.P_c;
        %P_c_current = LOCAL_CONSTS.P_c;

        %span_current = [P_c_current, 2500];

        %[p_m_vec, coef_m_mat] = pmvec_coefmmat_of_eigenmodes(m, LOCAL_CONSTS, span);
        %[p_m_vec, coef_m_mat] = pmvec_coefmmat_of_eigenmodes_analytically(m, CONSTS, span);
        p_m_vec = func_p_m(m, LOCAL_CONSTS, span, false);

        for ni = 1:size(p_m_vec, 1)

            z_det_vec(ni, i) = p_m_vec(ni);

        end

    end

    for k = 1:size(omega_vec, 1)
        for kk = 1:size(p_m_vec, 1)
            if(z_det_vec(kk, k)==0)
                z_det_vec(kk, k)=NaN;
            end
        end
    end

    for ni = 1:size(p_m_vec, 1)

        if (plot_data)
            legend_vec = [legend_vec; sprintf('n=%d', ni)];

            plot(h_dispers, omega_vec./omega_H, z_det_vec(ni, :), '.-', ...
                'color', [0.3; 0.0; 1.0*(ni-1)/(size(p_m_vec, 1)-1)] ...
            );
        end
    end


    if (plot_data)
        plot(h_dispers, omega_vec./omega_H, P_vec, 'm--', ...
                        omega_vec./omega_H, P_c_vec, 'r--');
        legend(legend_vec, 'P(\omega)'); xlabel('\omega/\omega_H'); ylabel('p_{mn}')
    end

end